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Emergence of nonlinear behavior in the dynamics of ultracold bosons 
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We study the evolution of a system of interacting ultracold bosons, which presents nonlinear, 
chaotic, behaviors in the limit of very large number of particles. Using the spectral entropy as an 
indicator of chaos and three different numerical approaches : Exact diagonalization, truncated Hu- 
simi method and mean-field (Gross-Pitaevskii) approximation, we put into evidence the destructive 
impact of quantum noise on the emergence of the nonlinear dynamics. 


PACS numbers: 03.75.Lm, 05.45.a, 03.75.Kk 

I. INTRODUCTION 

Ultracold atoms are clean, controllable, highly flexible 
experimental systems, that can often be modeled from 
first principles. For such reasons they have become in 
recent years a preferred testing ground for quantum 
many-body effects [1]. A particularly interesting example 
is the emergence of chaotic behaviors in quantum sys¬ 
tems : The Schrodinger equation is linear and hence can¬ 
not display chaotic behavior in the classical sense (i.e. 
chaos associated to a sensitivity to initial conditions), ho¬ 
wever, the classical world very often display chaos. Bose- 
Einstein condensates can be produced in laboratories at 
mesoscopic sizes (up to 10 8 atoms) [2] in intrinsically 
quantum-coherent states, they are thus ideal systems for 
the study of the transition form quantum to classical be¬ 
havior, and, in particular, the emergence of chaotic beha¬ 
viors. The corresponding quantum many-body problem 
(with binary contact interactions) can be treated by de¬ 
composing the matter-wave field ip in a macroscopic part 
describing a single particle wave function ip(x,t), called 
“condensed fraction”, plus a fluctuating quantum field 
corresponding to the excitations of the matter-wave field. 
The wave function ip(x,t) obeys the well-known Gross- 
Pitaevskii equation (GPE), which includes a nonlinear 
term [3-5] : 

= i Hl + 9 N l^(M)| 2 ) ip{x,t), (1) 

where g is the nonlinear parameter proportional to the s- 
wave scattering length, N is the number of particles and 
Hi is the one-particle Hamiltonian. The nonlinear term 
models the “mean effect” of particle-particle interactions, 
and is simply proportional to the spatial probability den¬ 
sity |^(a;, t)| 2 [26]. Replacing the many-body operator ip 
by a single-particle wave function ip which leads to the 
effective equation (1) implies neglecting quantum fluctua¬ 
tion, but preserves the long-range order supposed to be 
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an essential property of condensates. It turns out that in 
a large variety of interesting situations this rather rough 
approach gives a very good description of the condensate 
dynamics [3] . It can nevertheless be surprising that a li¬ 
near exact problem can be accurately modeled by a nonli¬ 
near effective equation which may present a qualitatively 
different dynamics (e.g. sensitivity to initial conditions 
and chaos). The mean-field approximation hides the “mi¬ 
croscopic” origin of the nonlinearity, and one is left with 
a situation similar to that encountered in statistical phy¬ 
sics : The macroscopic dynamics of a gas can present 
qualitatively different characteristics (e.g. irreversibility) 
from the microscopic dynamics of individual molecules 
(which is reversible). The averaging over the microscopic 
dynamics leads to a loss of symmetry and to a qualita¬ 
tively different behavior in the macroscopic scale. One 
can speculate that the nonlinear behavior of a conden¬ 
sate arises in an analog way : Averaging over the micro¬ 
scopic (quantum) dynamics produces a qualitatively dif¬ 
ferent behavior. A particularly spectacular manifestation 
of this is the “quasiclassical” chaos that can appear in the 
dynamics of a condensate [6-10], that is, chaos related to 
sensitivity to initial conditions. In the present work, we 
will use a simple model displaying quasiclassical chaos 
in the mean-field approximation and will compare it to 
the solutions of the many-body problem. By adequately 
choosing the representation of the dynamics we will show 
that one can observe “traces” of the quasiclassical chaos 
in the behavior of the many-body system even with a 
limited number of atoms. 


II. THE MODEL 

The toy model used here consists in a ultracold boson 
gas placed in an accelerated (or tilted) optical poten¬ 
tial [11, 12] corresponding to the single-particle Hamilto¬ 
nian 

1 d 2 

Hi = - ^— I r + Vocos(2nx) + Fx (2) 

7T 1 OX 1 

where we used normalized units [6] such that the lattice 
constant of the periodic part of the potential is 1, ener¬ 
gies are measured in units of the so-called recoil energy 
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( hoj r = H 2 k 2 L /2M : with k B = 7r /d where d is the lattice 
constant of the periodic part of the potential and M is 
the mass of the atoms) and we have set H = 1. In presence 
of interactions, a nonlinear term g\ip(x,t) | 2 is added to 
the above Hamiltonian leading to Eq. (1). We confine 
the dynamics to three adjacent sites s = —1,0,1 placed 
at x = —1,0,1 with Dirichlet boundary conditions [27] 
and restrict the dynamics to the lowest band of the sys¬ 
tem ; one can then write the solution of Eq. (2) in the 
form 

4>{x,t) = C-ie lFt w(x+l)+cow(x)+Cie~ lFt 'w(x— 1) (3) 

where w(x) is the Wannier-Stark function [eigenstate 
of (2)] of the fundamental ladder [28] associated to the 
x = 0 site [13]. In the case j / 0 , we keep the above 
form, but we let the coefficients c s depend on time. Then, 
inserting Eq. (3) in Eq. (1) one obtains a set of coupled 
differential equations that can be integrated numerically : 

ic s = Fsc s + UN\c s \ 2 c s 

+ J±N (2|c s | 2 c s ±i + |c s ^=i| 2 c s ^=i + c* ±1 c^) (4) 

with c _2 = C 2 = 0. The first term is simply the energy of 
site s, the second term [with U = g f dxw i {x)\ accounts 
for interactions of particles at the same site, and the two 
last terms correspond to the exchange of particles bet¬ 
ween neighbor sites, with [29] 

J+=g J dxw 3 (x)w(x + 1) 

J- =g J dxw(x)w 3 (x + 1). (5) 

It is useful to write the wave function components as 
amplitude-phase variables, defined by c s = \JT s e™ B . 

The dynamics of such system has been studied in pre¬ 
vious works [6, 8, 14, 15], where it was in particular sho¬ 
wed that quasiclassical chaos in the above system has the 
structure prescribed by the Kolmogorov-Arnold-Moser 
theorem, which becomes apparent in a Poincare section 
of the dynamics, Fig. la. We represented in Fig. lb the 
same dynamics described in terms of the spectral entropy 
S. Given a dynamical function f(t G [0,f max ]) (e.g. the 
average position (x) ( t )), the spectral entropy is related 
to the power spectrum F(i/) of f(t) by 

S—iiEbWlnFM, 

V 

where n v = 1 + t max /St is the number of frequency com¬ 
ponents of f(t) (with St the corresponding resolution). 
Roughly speaking, the spectral entropy gives the num¬ 
ber of significant frequency components present in the 
spectrum, and is thus a good indicator of (quasi)classical 
chaos. The aim of the present work, is to understand 
how the (linear) dynamics of the exact many-body pro¬ 
blem can approach the (quasiclassical) chaotic behavior 
observed in Fig. 1. 


III. NUMERICAL APPROACHES 

The exact many-body problem can be described by 
the Bose-Hubbard Hamiltonian [4, 16, 17]. The single¬ 
particle Hamiltonian corresponds to Eq. (2) , whereas bi¬ 
nary atom-atom contact interactions arise from a term 
(g/ 2) f dxifil (x)ifit (x)i/>(x)i/>(x) where 4>{x) is the matter- 
wave field, which is expanded in the Wannier-Stark basis 

= E s w(x — s)a s where a s is the boson annihilation 
operator for the site s. This leads to a Bose-Hubbard 
Hamiltonian 

Hbh = J2 S [ Fsn s + f n s {n s - 1) 

+ (J+a\a s+1 n s + J_a|a s _in s + h.c.)] 1 

where n s = aja s is the number operator for site s. Note 
that the two “kinetic energy” terms on the second line 
are different from those of the usual Bose-Hubbard Ha¬ 
miltonian as, contrary to Wannier functions of a perio¬ 
dic potential, Wannier-Stark functions are eigenstates of 
the one-particle Hamiltonian Eq. (2). Transport here is 
due to interactions between neighbor sites, not to tun¬ 
nel effect, and is thus proportional to the interaction 
strength g [cf. Eq. (5)]. One can convince oneself that 
the Bose-Hubbard Hamiltonian describes the many-body 
Wannier-Stark problem by noticing that the mean-field 
approximation which consists in replacing the operators 
a s by c-numbers c s leads to a classical Hamiltonian whose 
equations of motions correspond to the GPE, Eq.(4). In 
the following, we will compare the solution of the many- 
body Schrodinger equation 

= H B „ I*,*), (7) 

to the solution of Eq. (4). It is challenging to directly dia¬ 
gonalize Eq. (6) which is of dimension (IV+l)(A-|-2)/2 ~ 
TV 2 for atom numbers N larger than a few tens, so we 
used the so-called “Lanczos diagonalization” (LD) [18- 
20], which consists in using a truncated many-body evo¬ 
lution operator over a short time interval dt : 

n—0 

In the Lanczos scheme, in order to evolve the 
wave function \ipBH(t)} during dt one first builds 
an orthonormal basis |Zi),.., \l nK ) of the sub¬ 
space spanned by hk “Krylov vectors” defined as 
H 2 \ij) BH (Z)), .., H nK \ijjBn(t))- 
The resulting tridiagonal Hamiltonian of dimension 
tik ~ 10 — 20 <C N 2 can be then diagonalized to obtain 
the Lanczos eigenvectors |zq),.., \v nK ) and eigenvalues 
aq,.., a nK . The propagation of the wave function is then 
approximated as : 

n K 

I ipBH{t + dt)) = ^(i/ i |^ i3ff (t))e _l “ idt |r'i)- 

i= 1 
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Figure 1: Quasiclassical chaos in a tilted lattice, obtained by integrating the Gross-Pitaevskii equation (4). (a) Poincare 
section corresponding to 7_i = 0.1, do — 9 -1 = 0. (b) Spectral entropy calculated between t = 0 and t = 800. Parameters are 
V 0 = 5, F = 0.25, gN = 0.2. 


As the error in this approximation can be controlled 
by adjusting the time step dt, the Lanczos procedure is 
considered as an “exact diagonalization” technique. Each 
time step implies the diagonalization of a uk x hk ma¬ 
trix thus the method is only interesting if % <C TV 2 . The 
Lanczos procedure allows us to treat the many-body pro¬ 
blem with atom numbers TV < 500 in reasonable compu¬ 
ter times. 

Describing the system by a set of complex numbers 
c s instead of quantum operators a s , the GPE potentially 
neglects two important effects : (i) quantum noise, i.e the 
quantum fluctuations around the expectation value of a s 
(ii) quantum interferences, the fact that two different tra¬ 
jectories can interfere [30]. The truncated Husimi method 
(THM) allows one to reinsert quantum noise in a GPE 
description by propagating independently sets of initial 
conditions whose distribution reflects the quantum noise. 
However, the THM neglects quantum interference effects 
that may occur between the different trajectories and is 
thus an intermediary approach between the GPE and the 
exact solution of the many-body system particularly sui¬ 
ted to identify the role of the quantum noise [10, 21, 22). 
One thus needs to choose a rule relating the complex va¬ 
riables {c'J serving as initials conditions to the GPE to 
the many-body problem, i.e to choose a family of many- 
body states able to connect a many-body state |H ({cj)) 
and the corresponding set of GPE initial conditions {cj. 
With M sites (M = 3 in this work), one can choose 
the so-called SU(M) coherent states [23] with finite to¬ 
tal number of atoms TV, which are defined, in the Fock 
basis, as 

r n s 

|H({c s })) = v / M £ nf 

n 1 +..+n M =N s v n S' 

with J2 S l c s 1 2 = 1- The advantage of SU{M ) states, as 
compared to the usual Glauber coherent states, is that 

the total number of atoms TV is fixed : {( TV — (JV) ) ) = 


0, with TV = which is a situation closer to the 

that encountered in a real experiment. We then need 
to choose a quantity to characterize the quantum noise 
in the system. In this work we use the Husimi func¬ 
tion Qn which, given a coherent state |fi ({cj)), is sim¬ 
ply defined as its projection over all possible states : 
Qn ({c'J) = | (H' ({c'J)| H ({c s })> | 2 . The width of this 
function decreases as TV -1 / 2 and thus gives a good repre¬ 
sentation of the quantum noise. For a given initial many- 
body state |H{c s }), we generate the corresponding dis¬ 
tribution of one-particle initial conditions D (c s ) = {cj 
which mimics the Husimi representation of the |0{c s }). 
Each of these states is then propagated independently 
according to the GPE and observables are obtained by 
averaging over the distribution of final states. For large 
enough TV, the width of the Husimi distribution tends to 
0 and one recovers the results of the GPE. Hence, the 
THM allows one to use the GPE to simulate the evolu¬ 
tion of not so large number of atoms (including quantum 
noise but not interference effects), and it turns out to be 
a convenient tool to study of the emergence of chaos in 
our system as the number of atoms increases. 


IV. EMERGENCE OF QUASICLASSICAL 
CHAOS 

We compare in Fig. 2 the time evolution of the system 
- represented here by the average position of the wave 
packet - calculated according to the different methods 
described in sec. Ill, for a small TV = 30 [plot (a)] and 
a larger (but not large) TV = 500 [plot (b)[ number of 
atoms. In order to have a pertinent comparison with the 
GPE, gN is the same for all simulations. The GPE cal¬ 
culation displays a clear irregular behavior associated to 
the presence of quasiclassical chaos. No such behavior is 
observed for TV = 30 atoms either with the LD or the 
THM. In plot (b), with TV = 500, one sees that the GPE 
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Figure 2: Evolution of the average position of a wave packet 
calculated with different methods : Gross-Pitaevskii equation 
(cyan), Lanczos diagonalization (green), and the truncated 
Husimi method (red), for (a) IV = 30 and (b) N = 500. In 
panel (a), we also represented by blue circles the result obtai¬ 
ned by direct diagonalization of the full many-body problem, 
which is in perfect agreement with the result of the Lanc¬ 
zos method, illustrating the accuracy of the later. The initial 
state is a £1/(3) coherent state with c_i = \/0.5, Co = \AL25 
and ci = i\J 0.25. Other parameters are gN = 0.2, Vo = 7, 
F = 0.1. 


result “sticks” to the many-body calculations for times 
t < 150, and that the many-body calculations even dis¬ 
play the first oscillations observed in the GPE evolution. 

In order to observe the emergence of quasiclassical 
chaos as the number of particles increases, we calculated 
the spectral entropy for a complete set of initial SU( 3) 
states using both LD and THM. The width of the Hu¬ 
simi function corresponding to a given initial condition 
is determined, and we assume that for any initial condi¬ 
tion {c s }, the Husimi function can be approximated by 
a square distribution with same width centered around 
{c s }. We checked numerically that this approximation 
has no impact on our results. Figure 3 compares the dis¬ 
tribution of the spectral entropy computed by LD and 
by the THM for N = 30 and N = 400 atoms, and shows 
that not only the results become more similar when the 
atom number increases, but also that the characteristic 
features observed in the quasiclassical Poincare section 
Fig. la tend to emerge (this trend is confirmed for inter¬ 
mediate atom numbers). These results show that THM 
provides a good representation of the exact many-body 
behavior for large atom numbers and can advantageously 
used for comparisons with the quasiclassical behavior. 

In figure 4 we present a calculation of the variance 
of the wave packet position - which is an indication of 
the erratic character of the dynamics, that is quasiclas¬ 
sical chaos - by the three methods and for N = 30 (top 
row) and N = 400 (bottom row). The GPE picture (first 
column) is obviously independent of the atom number 
(provided gN is constant). Both LD (center column) and 



Figure 3: Spectral entropy distribution in phase space ob¬ 
tained by LD (left column) and by THM (right column) with 
N = 30 (top row) N = 400 (bottom row) and same parame¬ 
ters as Fig. 1. When N increases, the THM provides a quite 
good representation of the mean-field behavior for large N, 
and the phase-space structure of the many-body problem be¬ 
comes similar to the Poincare section obtained from the GPE, 
plot (a) of Fig. 1. 


THM (right column) instead show a clear dependence in 
N and are rather different, both in shape as in the am¬ 
plitude, from the GPE result, which shows that quasi¬ 
classical chaos is not fully developed for such atom num¬ 
bers. However, one can see however that the higher the 
number of atoms the closer the result is from the one 
obtained with GPE. The amplitude of the variances in¬ 
creases by a factor 2 for THM and by a factor 3.8 for 
the LD when the number of atoms is changed from 30 to 
400, and the surface of the chaotic zones (in yellow and 
red) also increases significantly. Thus, even if the atom 
numbers considered here are clearly too small to allow 
a definite conclusion, one can reasonably expect a much 
better convergence with the GPE for larger values of N. 


V. CONCLUSION 

In conclusion, we showed that the dynamics exhibited 
by the many-body problem as described by both exact 
diagonalization and the truncated Husimi method tends 
to converge to the quasiclassical chaos displayed by the 
mean-field approximation as the atom number increases. 
Our result suggests the existence of a “nonlinear characte¬ 
ristic time” increasing with some (monotonous) function 
of TV, during which the prediction of the mean-field ap- 
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Figure 4: Variance of the wave packet position calculated between t — 0 and t = 800 from the GPE (left column), LD (center 
column), and THM (right column) with N = 30 (top row) and N = 400 (bottom row) and the same parameters as Fig. 1. 


proximation and of the exact many-body problem agree. 
One may conjecture that this time is related to the ty¬ 
pical time for particle number fluctuations to affect si¬ 
gnificantly the dynamics, but we are presently unable 
to give a verifiable estimation of such time. In particu¬ 
lar, THM appears to approach very well the exact dyna¬ 
mics, which suggests that the essential effect leading to 
the emergence of quasiclassical chaos as the number of 
particles increases is the reduction of quantum noise. It 
would be interesting to analyze the influence of quantum 
interferences and diffusion as done e.g. by Carvalho et 
al. [24], who compared the phase space structure of the 
classical chaotic one-body delta-kicked harmonic oscilla¬ 
tor to the Wigner representation of its quantum coun¬ 
terpart, and showed that both decoherence and classical 
diffusion lead to a similarity of the two representations. 
Weiss and Teichmann [25], using a many-body system 
(with N = 1000) slightly different from ours, attribu¬ 
ted the observed differences between the mean-field and 
the many-body dynamics to the existence of entangle¬ 
ment in the many-body problem, and showed that these 
differences are reduced if decoherence is added to the 


system. The present work sheds a different light on the 
problem of the quantum-classical transition in a system 
displaying quasiclassical chaos, suggesting the quantum 
noise is the main cause of the differences. In order to 
reach more definitive conclusions, however, simulations 
of the many-body problem with atom numbers larger by 
(at least) one order of magnitude would be necessary. 
This is a fascinating yet challenging goal for future work. 
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